***This SAS program first generates the data for the nonlinear SEM and then does the fitting of it using PROC NLMIXED.
(The code for generating the data is more complex than necessary but can be used generally to simulate many datasets for a simulation study)
    

data param;
    reliab = .75;
	gam0= 0 ; gam1= 1 ; gam2= .5 ; gam3= 1; gam4 = .6; zeta1=.50 ;zeta2=.25; 

          b10 = 0 ; b20 = 0 ; b30 = 0 ; b40 = 0 ; b50 = 0 ; b60 = 0 ;
              b11 = .7 ; b21 = .3 ;
              b32 = .4 ; b42 = .8 ;
              b53 = .5 ; b63 = .4 ;

              mf1 = 0; 
              ph1= 1 ; 

       ps1 = (b11*b11*ph1*(1-reliab)/reliab);
       ps2 = (b21*b21*ph1*(1-reliab)/reliab);
       ps3 = (b32*b32*1.5*(1-reliab)/reliab);
       ps4 = (b42*b42*1.5*(1-reliab)/reliab);
       ps7 = (ph1*(1-reliab)/reliab);
       ps8 = (1.5*(1-reliab)/reliab);

       ps5 = (b53*b53*1.5*(1-reliab)/reliab);
       ps6 = (b63*b63*1.5*(1-reliab)/reliab);
       ps9 = (1.5*(1-reliab)/reliab);

     run;
proc print data = param; run;

%macro dosimulate(sam,obs);
     data start;
      do i=1 to %eval(&sam*&obs);
         %do j=1 %to 12;
         x&j=rannor(125);
         %end;
         output;
      end;

     data a;  set start;
     if _n_=1 then set param;

     f1=  mf1 + sqrt(ph1)*x1;
     f3=  gam0 + gam1*f1 + sqrt(zeta1)*x2;
     f2 = gam2 + gam3*exp(gam4*f3) + sqrt(zeta2)*x3; 

	 ****f2 is eta1;
	 ****f3 is eta2;


     Y1   = b10+ b11*f1   + sqrt(ps1)*x4;
     Y2   = b20+ b21*f1   + sqrt(ps2)*x5;
     Y3   = b30+ b32*f2   + sqrt(ps3)*x6;
     Y4   = b40+ b42*f2   + sqrt(ps4)*x7;
     Y5   = b50+ b53*f3   + sqrt(ps5)*x8;
     Y6   = b60+ b63*f3   + sqrt(ps6)*x9;
     Y7   =          f1   + sqrt(ps7)*x10;
     Y8   =          f2   + sqrt(ps8)*x11;
     Y9   =          f3   + sqrt(ps9)*x12;
output;
     keep Y1 Y2 Y3 Y4 Y5 Y6 Y7 Y8 Y9 f1 f2 f3;
     run;

     data a; set a;
          %do p=1 %to &sam;
          if %eval(1+(&p-1)*&obs)<=_n_<=%eval(&p*&obs) then iter=&p;
		  count = _n_;
		  dummy = 1;
          %end;
		  run;

          proc means data=a noprint; by iter;
          var y1-y9;
          output out=outmean mean=mean1-mean9;
          run;

          data a; set a;
           merge a outmean; by iter;
          run;

%mend dosimulate;

%dosimulate(1,500);


/*data forwinbugs; set a; file 
"C:\file\verynonlin.dat";
put y1-y9;
run;*/

*************************************Using nlmixed to fit the data***************;

proc nlmixed data = a (keep = dummy count y1-y9 iter where = (iter=1)); by iter;

***starting values;

****here I am using .5 * variance of the observed variable as the starting value of of psi, everything
else gets (by default) a starting value of 1;
parms
psi1 = .30 psi2 = .06 psi3 = .129 psi4=.53
psi5 = .24 psi6 = .15 psi7 = .64 psi8 = .80 psi9 = .96

bounds psi1-psi9>=0, phi1>=0, ddelta1-ddelta2>=0;
 
 f3i = gam0 + gam1*f1i + delta1i;
 f2i = gam2 + gam3*exp(gam4*f3i + gam5*f1i) + delta2i; 

     mu1   = b10+ b11*f1i   ;
     mu2   = b20+ b21*f1i   ;
     mu3   = b30+ b32*f2i   ;
     mu4   = b40+ b42*f2i   ;
     mu5   = b50+ b53*f3i   ;
     mu6   = b60+ b63*f3i   ;
     mu7   =          f1i   ;
     mu8   =          f2i   ;
     mu9   =          f3i   ;

ll_factpart = -.5*log(psi1) - (1/(2*psi1)) * (y1 - mu1)**2
			  -.5*log(psi2) - (1/(2*psi2)) * (y2 - mu2)**2
			  -.5*log(psi3) - (1/(2*psi3)) * (y3 - mu3)**2
			  -.5*log(psi4) - (1/(2*psi4)) * (y4 - mu4)**2
			  -.5*log(psi5) - (1/(2*psi5)) * (y5 - mu5)**2
			  -.5*log(psi6) - (1/(2*psi6)) * (y6 - mu6)**2
			  -.5*log(psi7) - (1/(2*psi7)) * (y7 - mu7)**2
			  -.5*log(psi8) - (1/(2*psi8)) * (y8 - mu8)**2
			  -.5*log(psi9) - (1/(2*psi9)) * (y9 - mu9)**2;

model dummy ~ general(ll_factpart);

random f1i delta1i delta2i ~ normal([mf1, 0, 0], [phi1, 
                                      0, ddelta1, 
                                      0,     0,   ddelta2]) subject = count;

run;